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A METHOD FOR ANALYZING DYNAMIC STALL 


OF HELICOPTER ROTOR BLADES 

By Peter Crimi and Barry L. Reeves 
Avco Systems Division 


SUMMARY 


A method has been developed for analyzing the dynamic stall of helicopter 
rotor blades. The method employs a model for each of the basic flow elements 
involved in the unsteady stall of a two-dimensional airfoil in incompressible 
flow. The interaction of these elements is analyzed using a digital computer. 

Calculations of the loading on an airfoil during transient and sinusoidal 
pitching motions are in good qualitative agreement with measured loads. Dynamic 
overshoot, or lift in excess of the maximum static value, and unstable moment 
variation are in clear evidence in the computed results. Quantitative differ- 
ences can be attributed in part to the use of a linearized representation of 
the potential flow and a quasi-steady model of the viscous mixing region. 

Computations were also performed of the loading and pitch response which 
result from unsteady stall induced by a series of discrete vortices convected 
past an elastically restrained airfoil. The results were used to confirm that 
large torsional response of helicopter blades during a maneuver which had been 
detected in flight tests can be attributed to dynamic stall induced by pre- 
viously formed tip vortices. 




INTRODUCTION 


Helicopter operation at high forward speed requires that the retreating 
blade develop a high lift coefficient because of the diminished dynamic pres- 
sure over that part of the rotor disc. As a result, under conditions of maxi- 
mum performance, the flow periodically separates from and reattaches to each 
blade, giving rise to severe oscillatory control loads, increased vibration 
levels and, at times, a torsional aeroelastic instability. Consequently, un- 
steady blade stall seriously limits helicopter performance (Refs. 1 and 2). 

Data from wind tunnel tests of oscillating airfoils have shovm that stall- 
ing under unsteady conditions differs markedly from static airfoil stall (Ref. 
3). Lift coefficients 20 to 30 percent in excess of the maximum static value 
have been measured. Also, there is a substantial loading hysteresis accompanied 
by large nose-down aerodynamic moments varying in time in such a way as to ex- 
tract energy from the free stream. 

There have been a number of analytical studies of the problem in recent 
years (Refs. 4, 5 and 6). The results are of limited applicability, however, 
because of the use of empirical methods. 

The study reported here was directed toward developing a general method 
for analyzing dynamic stall and applying the method in the analysis of a type 
of wake- induced stall detected in helicopter flight- test data. The method 
developed utilizes mathematical representations of each of the flow elements 
involved in unsteady airfoil stall. The interactions of these elements are 
analyzed using a digital computer. Consideration has been limited to two- 
dimensional, incompressible flow. However, the method has been formulated to 
account for arbitrary prescribed airfoil motions and airfoil section charac- 
teristics. The method is, therefore, applicable to a number of aerodynamic 
and aeroelastic problems involving unsteady airfoil stall, such as stall 
flutter of propellers, rotors and compressor blading, rotating stall in axial- 
flow compressors and response of aircraft to severe gusts. 
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SYMBOLS 


a mn » ^mn » 
a kn» ^kn 

B n 

b 

C 


n 

H 

H 

h 


J 

k 

k 0 

1 

1 

s 

M 


m 


coefficients in series representation of y 
distance aft of midchord of pitch axis in semichords 
coefficients in linear equations for A n ’s and B n ’s 

coefficients in representation of a 
airfoil semichord, m 
section camber distribution, m 
lift coefficient, = l/(pU 2 b) 
moment coefficient, C m = M/(2 p U 2 b 2 ) 
normal-force coefficient, C n =N/(pU 2 b) 
pressure coefficient, C p = 2 (p-p^ ) /( pU z ) 
airfoil chord, m 

coefficient in series representation of C 
shear layer shape factor, H= 0/8 * 
vertical separation of vortices, in semichords 
section plunging velocity, m/s 
shear layer parameter, J = 6*/ 8* 

reciprocal turbulent Reynolds number, = * m /u e 6 
reduced frequency, k = cob/U 

dimensionless pitch natural frequency, = <*q b/U 
lift per unit span, N/m 
length of dead-air region, m 
Mach number 

moment per unit span about pitch axis, positive to increase 0 p 
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Nf 


N 

V 

N a 

N 



r m> ? k 



T 

TM 


t 

t 

t 

n 

U 


number of terms in series representations of T and C 
number of wake elements 

number of coefficients in series representation of y 
number of grid points in the 7/ -direction 
number of source elements 

force component normal to chord line per unit span, N/m 
pressure, N/m or psi 

2 

dimensionless pressure, p = p/p U Q 

magnitude of flow external to boundary layer, m/s 

ratio of streamwise boundary-layer flow component to U Q 

ratio of normal boundary-layer flow component to U Q /\/Rej^ 

Reynolds number based on length indicated by subscript 

integral over shear layer, 

R = (2 U 0 /Kfl H u e ) J ( r / p u e 2 ) ( d q s / dy ) dy 

0 

inhomogeneous terms of equations for A n ’s and B n ’s 
leading-edge radius , m 

streamwise vortex spacing, in semichords 
streamwise boundary-layer coordinate 
m th streamwise boundary-layer grid point 
section thickness distribution, m 
rotor blade torsional moment, lb -in 
time, s 

dimensionless time, t = U Q c/b 
coefficient in series representation of T 
instantaneous free-stream speed, m/s 
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* 

u 

V 


w 

<x, y) 


(x l • n } 




Y 


a 


y 


y\j 

} % 

ACp 

Ap 

Ac 

A 0 
8 


reference velocity, m/s 

flow perturbation in x-direction, m/s 

flow external to viscous mixing region 

root-mean-square fluctuation of free-stream due to turbulence 
flow perturbation in y -direction, m/s 

downwash due to airfoil incidence, motions and camber, m/s 

coordinate system with origirf at midchord of mean position 
of airfoil 

coordinate system with origin at leading edge of mean position 
of airfoil 

value of x at the terminus of the vortex wake 
coordinate of n th vortex wake element 

coordinate of n th point at which flow-tangency condition is 
imposed 

coordinate of i th source element 

location of midpoint of i th source element 

ordinate of airfoil surface with respect to chordline 

angle of attack, rad or deg 

bound vortex strength, . m/s 

wake vortex strength, m/s 

wake vortex strength at n tJl wake point, m/s 

CP L “ C PU 
PL “Pu 

increment for integration in time, s 

half amplitude of sinusoidal pitching, rad or deg 

boundary layer thickness defined as distance from surface 
at which u o q s = 0.995 q e , m 
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boundary layer or shear layer displacement thickness, m 


SVRi^/b 
S* VRT b /b 

eddy viscosity, m 2 /s 

dimensionless boundary laye r coordinate, ratio of distance 
from the surface to b/ V Re~ 

n th normal boundary-layer grid point 

boundary layer or shear layer displacement thickness, m 
pitch angle, rad or deg 

mean pitch angle or pitch angle of zero restoring moment, deg 
shear layer energy thickness 

pressure gradient parameter, A = - ( 8 2 dp/dx)/(nq e ) 
viscosity, N-s/m 2 
kinematic viscosity, m 2 /s 
1 + e/ v 
density, kg/m^ 

strength of source distribution representing dead-air region, 
m/s 

total source strength representing airfoil thickness and dead- 
air region, m/s 

2 

shear stress, N/m 

2 

perturbation velocity potential, m /s 

rotor blade azimuth angle, measured from downwind direction, 
deg 

angular frequency of oscillation, rad/s 
pitch natural frequency, rad/s 


Subscripts: 

b 

d 

L 

lam 

R 


s 

c 

turb 


U 


oo 


point at which pressure recovery begins in viscous mixing 
region 

dead-air region 

lower surface of airfoil 

laminar 

point of reattachment of turbulent shear layer 
separation point 
transition point 
turbulent 

upper surface of airfoil 
free stream 




STALL MECHANISMS AND FLOW ELEMENTS 


The mechanisms of stall onset are extremely complex and depend on many 
parameters, including Reynolds number, Mach number, leading-edge radius, thick- 
ness, camber, sweep, and the pressures imposed by unsteady motion. It is 
generally accepted that a given airfoil in steady flow stalls in one of three 
ways (Ref. 7). The three types of stall are trailing-edge stall, leading-edge 
stall and thin airfoil stall. Examples of each of these types are documented, 
with measured pressure distributions and boundary-layer profiles, and discussed 
in detail in Ref. 7. 

Trailing-edge stall is the most easily identified of the three types, 
being due to the separation of the turbulent boundary layer near the trailing 
edge. Increasing incidence moves the point of separation progessively forward 
along the airfoil, resulting in a gradual decrease in lift and increase in 
drag, as indicated in Figure la. This type of stall generally occurs on rela- 
tively thick airfoils at high Reynolds numbers. 

Leading-edge stall is related to the formation of a small separation bubble 
near the leading edge. At a fairly low incidence, laminar separation occurs 
near the point of minimum pressure at the leading edge. The flow reattaches a 
short distance downstream of the separation point because of transition from 
laminar to turbulent flow in the free shear layer with subsequent turbulent 
mixing and reattachment. As the angle of attack increases, the bubble moves 
closer to the leading edge, grows slightly shorter and somewhat thicker. The 
bubble has almost no effect on integrated loads, because it is never more than 
a few percent of the chord in length. At some angle of attack, the bubble 
bursts and the flow separates from the entire upper surface of the airfoil, 
resulting in a sudden loss in lift, as indicated in Figure lb. The precise 
reason for the bursting of the laminar bubble has been the subject of consider- 
able controversy. There have been correlations attempted with bubble length 
and with boundary-layer momentum thickness at the point of laminar separation, 
with little success. There is a strong indication, though, that there is some 
maximum amount of pressure recovery which can occur in .the turbulent mixing 
zone and still allow reattachment, and at some incidence the required recovery 
exceeds this maximum, causing sudden 'separation. A thorough and well-ordered 
discussion of the various theories and evidence related to leading-edge stall 
is given in Ref. 8. 

Thin-airfoil stall, which occurs at relatively low Reynolds numbers on 
thin airfoils, is characterized by the appearance of a laminar bubble spring- 
ing from the leading edge at a relatively low incidence. Unlike the bubble 
formed prior to leading-edge stall, its point of separation remains fixed with 
increasing incidence while the bubble grows progressively larger. The processes 
of bubble formation and reattachment are not well understood (Ref. 8). The re- 
sulting lift curve is as sketched in Figure lc. Because of the uncertainty as 
to the precise mechanism of thin-airfoil stall and its relatively infrequent 
occurrence, it was decided not to attempt to model the thin-airfoil stall 
mechanism. Elements necessary to account for both leading-edge and trailing- 
edge stall have been represented in the method, however. 
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Angie of attack 


a Trailing - edge stal I 




b Leading -edge stall 



c Thin- airfoil stall 


Figure 1 THE THREE TYPES OF AIRFOIL STALL 
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The primary flow elements of unsteady leading-edge or trailing-edge stall 
of a two-dimensional airfoil can be specifically identified, as indicated in 
Figure 2. When the flow is attached (Figure 2a), the flow elements are: 

(1) a laminar boundary layer extending from the stagnation point over the 
leading edge; (2) a leading-edge separation bubble (if separation occurs prior 
to transition); (3) a turbulent boundary layer from the reattachment point of 
the leading-edge bubble (or the transition point) to the trailing edge; and 

(4) a potential flow over the airfoil, including the effects of a vortical 
wake generated by the variation in time of the circulation about the airfoil. 

If the airfoil undergoes leading-edge stall (Figure 2b), the flow elements 
are: (1) a laminar boundary layer to the point of separation; (2) a laminar 

constant-pressure shear layer to the point of transition; (3) a turbulent 
constant-pressure shear layer; (4) a turbulent pressure-recovery region; and 

(5) a potential flow over the airfoil and external to the viscous mixing 

region, again including a vortical wake. If trailing-edge stall occurs (Fig- 
ure 2c), the flow elements are: (1) the laminar boundary layer; (2) the lead- 

ing-edge bubble (if laminar separation occurs prior to transition); (3) the 
turbulent boundary layer; (4) a turbulent constant-pressure shear layer; (5) a 
turbulent pressure-recovery region; and (6) a potential flow with vortical wake. 

It was felt that, with the analytical tools at hand, a practical method 
could not be developed which precisely modeled each of these flow elements. 
Therefore, quantitative capability was sacrificed where approximations could 
be made while still retaining essential qualitative features of the flow. 
Specifically, approximations were employed in modeling the potential flow, the 
viscous mixing regions and the leading-edge bubble. The potential flow is 
derived from linearized boundary conditions. The viscous mixing regions and 
the leading-edge bubble are analyzed assuming quasi-steady flow. Also, inter- 
actions of the dead-air region with the inviscid flow are only indirectly taken 
into account. On the other hand, it could not be determined a priori what 
approximations would be permissible in computing the flow in the boundary 
layers, so a complete second-order, unsteady finite-difference method has been 
developed to analyze the boundary layers. The formulations used to represent 
the flow elements are given in the next section. 
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1. Laminar boundary layer y Laminar boundary layer 1. Laminar boundary layer 

2. Leading-edge bubble 2 - Laminar mixing region 2. Leading-edge bubble 

3. Turbulent boundary layer 3 - Turbulent mixing region 3. Turbulent boundary layer 

I 4. Airfoil and vortex wake 4. Turbulent reattachment region 4. Turbulent mixing region 

, 5. Airfoil and vortex wake 5. Turbulent reattachment region 

6. Airfoil and vortex wake 

\ a> Attached flow b * Leading-edge stall c. Trailing-edge stall 


Figure 2 FLOW ELEMENTS 


REPRESENTATIONS OF FLOW ELEMENTS 


Potential Flow 

Given the airfoil section characteristics and motions, together with the 
distribution of pressure in the dead-air region if the airfoil is stalled, the 
flow and pressure over the airfoil must be determined in order to compute the 
integrated load and analyze the boundary layer. The problem is formulated as 
follows. 

Consider an airfoil of infinite span and chord 2b subjected to a uniform, 
incompressible free stream of magnitude U(t), as indicated in Figure 3. Let 
0 (t ) and h (t ) denote the pitch angle and plunging rate of the airfoil, as 
shown in Figure 3. Further, let camber and thickness distributions C (x ) and 
T(x), respectively, be defined by 

C(x) = -L [Y v (x) + Y l (x) ] 

TOO = -L [ Y v 00 - Y l (x) ] 


where Y u and Yl are the ordinates of the upper and lower airfoil surfaces, 
respectively, measured from the chord line. 

Only the problem of a stalled airfoil need be considered, because the 
solution for attached flow is readily recovered as a special case. The co- 
ordinates of the separation and reattachment points and the prescribed pres- 
sure between those points are denoted x g , Xr and p^ , respectively. 

It is assumed at this point that perturbations to the flow caused by the 
airfoil are small compared to the free-stream speed, so that second-order quan- 
tities are negligible and boundary conditions can be imposed on the x-axis. 

This assumption, which forms the basis of classical thin-airfoil theory 
(Ref. 9), is clearly questionable at angles of attack sufficient to cause 
stall. However, the solution so obtained will still retain the essential 
qualitative features of the physical flow. That is, the computed pressure 
distribution will depend on the imposed boundary conditions, as derived from 
the foil motions, in approximately the same way as does the actual pressure 
distribution. At the same time, the solution of the linearized problem is 
much more readily incorporated in a large digital- computer program than would 
be the solution of the exact nonlinear problem. 

Thus, let u (x , y , t ) and v (x, y , t ) denote components of the flow pertur- 
bation in the x and y directions, respectively, and let p (x , y , t ) denote the 
pressure. The boundary conditions which must be satisfied, for x^ > b , are 
then given by: 
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± UT'- w , 
- UT'- vr , 


- b < x < x s ; 
x s < X < b ; 


+ 

v (s, o , t ) = 

v (x, 0~, t) = 


p (x, 0 + , t ) = p d (x, t ) , x s < x < b ; 

+ 

p (x, 0 > c ) = Pd (*, t ) , b < X < x R . 

where 

w = U - CO + h + (x - ab) 
and similarly for the case with x R < b. 


The (+) and (-) signs affixed to the zero arguments designate the sign of y in 
the limit y -* 0. It is further required that the flow be continuous at the 
trailing edge (the Kutta condition) and at the separation point. 


Assuming the flow is irrotational, a perturbation velocity potential <f> can 
be formulated by distributing source and vortex singularities on the x-axis: 


<f> - - 


1 

2 77 



(£, O tan 1 



- % ( £ O a/(x- 0 2 + y 2 


d£ 


where x Q marks the location of the starting vortex. It is convenient to 
separate out thickness effects. Thus, let 

a o (x, t) = 2UT / , -b < x < x s 

= 2UT + a (x, t) , x s < x < b ; 

= a (x , t ) , b < x < x R ; 

= 0 , x R < x < x Q . 

Also, the vortex strength for x>b is known in terms of the strength on the air- 
foil, through conservation of circulation: 


y (b, t) 


1 d 
u~ d7 



y (x, t ) dx 


( 1 ) 


and the vorticity downstream of the trailing edge is convected at the instan- 
taneous free-stream speed (to first order), retaining the strength it had upon 
formation at the trailing edge, as prescribed by Eq. (1). 
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With the potential formulated in this way, the two unknown functions are 

the vortex strength y on the interval (-b, b ) and the source strength o over 

the interval (xs»x R ). The boundary conditions, "written in terms of these 

two functions, are found to be (for x > b). 

R 


*£ 


cr (x, t ) + 


y(£,t) d£ 
X - f 


= 2 w — 




y w ( t ) d £ 
x- 7 


, - b < x < x s ; 


( 2 ) 


’ l b 


y(£,t)d£ l 

— = 2w 

X — £ 7T 


/ 


y w (£,e)d£ (3) 

, x s < x < b ; 


y (x, :) + — 

7T 


rf 


2fPoc-Pd> 2U / T'(£)d£ 


x- £ 


pU 


b 


X - £ 




+ 2 


- / 


cr(€,t) in |x - £ | d<f 


T'U) In |X- { I d*j , 


(4) 


x < x < b ; 
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1 


g(£, t)d£ 

x - £ 


2 (Poo-Pd> 

pu 


TXa££ 

s - ^ 


7T 



2 U 

77 



± JL)JL 

XJ dt ) 77 



0-( f , t) &i I x 


£|d<f 


+ 


U 

77 



T ' (a ea|x- ^ I d£ 


b < x < sj^- 


(5) 


The Cauchy principal value is taken for the singular integrals- A subscript 
w has been affixed to y for x > b to indicate that it is known in terms of y for 
x <b. 


Equations (2) through (5) were solved by first casting y in series form, 
with unknown coefficients, and taking a (with a singular term subtracted) to 
vary linearly between prescribed points on the x-axis, the value of a at each 
of those points being unknown . Equations (2) through (5) provided the rela- 
tions needed to solve for the unknown coefficients and the values of a . 

This procedure required that the functional form of y and o , i-e-, the 
locations and types of singularities, first be determined. For that purpose, 
an analytical solution was derived for the case of steady flow about a stalled 
flat plate with constant pressure in the dead-air region. The derivation is 
outlined in Appendix A. 

The solution for steady flow, as given in Appendix A, shows, first, that 
if is greater than b , y has the same functional form as for attached flow, 
i.e., y has a square-root singularity at the leading edge and is continuous 
at both x = x s and x = b. If x R is less than fc>> however, y has a square- root 
singularity on the downstream side of x = x R . Furthermore, regardless of the 
value of x R , a is zero at x = x s and has a square-root singularity at x = x R 
. (necessarily on the upstream side) . 
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Thus, y was taken to be of the form 



( 6 ) 


where cos 6 = x/b , U Q is a reference velocity and the coefficients A 0 through 
Atv^ and A R are functions of time. The term with the factor A R is dropped if 
either x R > b or x < x R . The factors involving x in the last two terms make 
Eq. (6) satisfy Eq. (1) identically, as can be verified by direct substitution. 


The form of a was defined by first dividing the interval (x s » x R ) into 
N a segments with end points x a> given by 


1 1 

x 05 = T (x R + x s> - T (x R " x s> cos 


(i - 1 ) if 
i = 1, 2, 


N ff + l; 


whereupon 


o (x, t ) 


U 


O 



(x ff . + l - x)B| + Csc-Xff.) B i + 1 


7 i +1 


(7) 


< X < X c 


i+1 


(Bl 


3 N 


o + 1 


0) 


where the B-*s are unknown functions of time. 


The coefficient A R in y , which is one of the unknowns when x R < b , is in 
fact proportional to B 0 • If Eqs. (6) and (7) are substituted in Eq. (3) and 
the limit x -> x R is taken, it is found that the result can only be finite at 
x = Xr (which it clearly must be) if 
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Thus, there are just N a 4- + 1 quantities to be determined at a given instant 

to completely define the flow and loading on the airfoil. 


A set of linear algebraic equations was derived, using Eqs. (2) through 
(5), to provide the relations needed to solve for the unknowns. Specifically, 
after substituting Eqs. (6) and (7) for y and a, respectively, the flow- 
tangency boundary condition (Eqs . (2) and (3)) was imposed at points x , where 

r- 7 n 


■y n 


b cos 


(n - 1) n 

n" 


n = 1, 2, • • • , Ny + 1 ; 


and the pressure boundary condition (Eqs. (4) and (5)) was imposed at points 
x , where 

V 

*•« ■ -T + *»« + l 


The resulting linear algebraic equations are of the form 


a mn A n + ^ b mn B n = r m , m = 1, 2, • . • , N y + 1 ; (8a) 

n =0 
n^l 

S kn A n + ^ b kn B n = ? k , k = 1, 2, • • • , N a (8b) 

n = 0 
n =/ 1 

Equations (8a) were derived from Eqs. (2) and (3), and Eqs. (8b) from Eqs. (4) 
and (5). The expressions obtained for the coefficients and inhomogeneous terms 
of Eqs. (8a) and (8b) are given in Appendix B. In the development of those 
relations, a number of assumptions and approximations were required. 

Derivatives with respect to time in Eqs. (2) through (5) were approximated 
by second-order finite differences. Given the value of, say, A q at times t, 
t - At and t - 2 At, the relation 

3 A 0 (t) -4A c (t - At) + A 0 (t - 2 At) 

A 0 (t) = — + O ((At) 2 ) 

° 2 At 

was utilized, and similarly for the other required time derivatives. 
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Thickness and camber distributions were represented by trigonometric series 
so the integrals in xdiich those functions appear could be evaluated analytically. 
Specifically, T and C were written in the form 


T 

T 

max 


Nf 

( 1 — cos 8 ) 

n = 1 



t fl sin n 8 


(9) 


C 

c 

max 


sin 6 



( 10 ) 


where cos 6 = x/b and T max and C max are the maximums of T and c, respectively. 
It was found that, using 24 terms, these series approximated the actual offsets 
of conventional NACA sections to within one percent (about .1 percent of chord) 
over the whole chord. The coefficients t Q and c n are computed from 


t 


n 


c 


n 


2 


77 T 


max 


2 


77 C 


max 



T ( 6 ) sin n 8 
1 — cos 8 


C ( 8 ) sin n 8 
sin 8 


d 8 


The wake dox^nwash integrals appearing on the right-hand side of Eqs. (2) 
and (3) x^ere evaluated by first specifying that the value of y w be known at 
discrete points along the x-axis, vrith the points spaced at a distance cor- 
responding to the time increment At used for integration in time. The strength 
at each point is obtained as follows. Numbering the wake points from 1 to N w 
(N w increases by one after each time step), x?ith point 1 at the trailing edge 
and point N w at x = x Q , then y w ^ (t) is computed from Eq. (1), using a second- 

order finite-difference approximation for the time derivative, while y w ^ (t) = 

y w (t-At) for n = 2, 3, ...,N . The vortex strength is assumed to vary 

n— 1 w 

linearly betx^een the designated points at xjhich y w is known. Note that the 
wake dox<mxjash integral contains a term involving unknown coefficients A (t) 
and A^(t) x^hich must be taken into account in the expressions for a mo and a^. 
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The solution of Eqs. (8a) and (8b), which is obtained by a simple elimina- 
tion scheme, provides the quantities needed to evaluate the pressure and flow 
at the airfoil surface. To first order, the flow at the surface is equal to 
U + u , and 


+ 

U (x, o , t) 


± 


y(x, 0 
2~ 


+ 


1 

2 77 



*(£ Odf 
x- € 



T^(£)d£ 
x - £ 


This expression clearly cannot be used to define the flow external to the 
boundary layer, because y is infinite at the leading edge. This difficulty 
was resolved by employing a correction factor derived by Lighthill (Ref. 10) 
which makes the first-order solution uniformly valid over the whole airfoil 
surface. Specifically, if q e denotes the flow at the surface, defined posi- 
tive in the clockwise sense as viewed in Figure 3, the Lighthill result gives, 
in the notation of this report, 


q e 0 = ± 


b + x 

b + x + r 0 /2 


1/2 


[U + u (x, o 


, t 


)] 


(ID 


where r Q is the leading-edge radius and the (+) and (-) signs apply to the 
upper and lower surfaces, respectively. Note that since r Q is of order 
q e differs from ±(U+u) by an amount which is of second order except in the 
immediate vicinity of the leading edge, and that q e is finite at the leading 
edge. The corresponding expression for the pressure coefficient, obtained 
from the linearized Bernoulli equation and Eq. (11), is 
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Explicit relations for q e and obtained by substituting Eqs. (6), (7) 
and (9) in Eqs. (11) and (12), are given in Appendix B. Similar expressions 
for lift and moment coefficients can be developed. However, those coefficients 
are more readily evaluated by numerically integrating Cp , so explicit expressions 
were not derived. 
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Boundary Layer 


Because the relative importance of the individual elements of the boundary- 
layer flow as they affect dynamic stall could not be established in advance, the 
representation of the boundary layer was made as general as possible to ensure 
that the effects of all essential elements were taken into account. The method 
of finite differences, rather than an integral method, was selected for the 
analysis of both laminar and turbulent boundary layers because accuracy and com- 
puter requirements are readily controlled and the simpler formulation allowed 
rapid initial check-out of computer coding. The Smith-Cebeci eddy-viscosity 
model (Ref. 11), which has given good results for a wide range of Reynolds 
numbers and various types of pressure distribution, was chosen to represent 
the turbulent shear. 

The boundary- layer equations were cast in the following dimensionless form 
for the cases of both laminar and turbulent flow: 

d % 9 % d< is 9 P 9 (_ 9 % \ (13) 

dt + qs <9s + = ~ ~ds~ + Tif V e TIT ) 

9 % 9 *r, < 14 > 

IT + — = 0 

Physical quantities relate to the variables in the above equations as follows. 
The flow component parallel to the surface is U Q q s , the flow component normal 
to the surface is U D /\/Re^, where Re b is Reynolds number based on semichord, 
time is bt/U Q , distance along the surface is bs , distance normal to the surface 
is b 17 / V Re^ and pressure is p u 2 p. For laminar flow, v e - 1 and for turbulent 
flow ’ ° ' 


1 


+ 


€ 

V 


where v is kinematic viscosity and c is eddy viscosity. In terms of the 
variables defined above (see Ref. 11) 


— = .16 7 7 2 Vr^ ( 1 - e~ K )‘ 

If ^ 
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and 


.0168 (qg/Uj 8 * 

VRi" b [ 1 + 5.5 (77/S) 6 ] 


= £ o ’ 1 t> Vo 


where 


5 * = 


f 


d7 7 


while d is the value of 77 at which q = .995 q /U and 77 is the value of 77 at 
, , , _ __ e o 1 o 1 

which e. = e . 

1 o 

The method employs variable step size in both the s and 77 directions. 
Finite-difference approximations for the derivatives appearing in Eq. ( 13 ) 
were derived from Taylor series. The error in each approximation is of the 
order of the square of the step size. The differences were formulated to 
allow computation of the flow at time t and streamwise coordinate s m4 q , given 
the flow at times t- At and t-2At at all streamwise coordinates and the flow 
at streamwise coordinates s m and for time t . The finite-difference approxi- 
mations used for dq s /dt , dq s /ds , dq s /d77 , d 2 q s /drj 2 and dp/ds are given in 
Appendix C. 

The solution for the flow at s ra + i is obtained by iteration. For the first 
iteration, the difference equations are linearized by employing the following 
approximate relations: 
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where d n , and f n are defined in Appendix C and 
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For the second and succeeding iterations, the values computed in the previous 
iteration are employed in place of the above extrapolations. 

Substitution of the finite- difference approximations in Eq. (13) then 
yields a set of linear algebraic equations of the form 


n ^ s m+l , n -1 + ^ s m + 1 , n ^m + l , n + 1 + ^ s m+l , n +2 


- F„ 


(15) 


n = 2, 3, • • • , — 1 


where the coefficients A n through D n and F n only involve flow quantitives com- 
puted at s m and at time t and at s m + j at previous time steps. After 

setting q p equal to zero and q o and q o equal to q^ /tJ , 

&4s m + l,l 4 4s m + l,N„ 4 s m+l, N„+l 4e m + l ° ’ 

Eqs. (15) are readily solved by successive elimination. Equation (14) can then 
be used to compute q , from 


q 


1 



drj 


using a trapezoidal approximation. 

The aforementioned iteration uses wall shear as the criterion for con- 
vergence, the allowable error being .1 percent. A maximum of five iterations 
is allowed at each s station. Some difficulty with convergence was encountered 
under extreme conditions. This problem was resolved by smoothing the eddy 
viscosity variation with 77 , as was done in Ref. 11. After 7 is computed, its 
value is replaced by the mean of the values at three adjacent points. 

The grid in the ?/ direction is computed from a geometric progression 

similar in form to the one used in Ref. 11. Given values of r and V2^l s 

■q is computed from ^ 

n 

7? n = (1 + r^) r ]a _ 1 - i? n _ 2 , n = 3, 4, • • • , + 1 

At a given instant, the boundary layer thickness at the trailing edge can be 
an order of magnitude larger than at the stagnation point. To accurately com- 
pute the flow at these extremes with a single grid in the r\ direction would 
require several hundred mesh points in that coordinate. About two hundred 
points are needed in the s-direction for a typical airfoil. Thus, an 
inordinate amount of computer storage would then be required, because the 
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value of q s at every mesh point for two time steps must be stored. Therefore, 

it was decided to make the rj - scale variable. If, upon completion of the flow 

computation, the boundary- layer thickness 8 exceeds % ^ is increased by 

7? — A 1 


a set amount, a new, expanded 77 -scale is computed and the flow quantities at 
the new mesh points are assigned by interpolation. The variation of the 
interpolated quantities is smoothed by three-point averaging, and computations 
proceed as before. The value of is stored with the flow quantities so the 
boundary- layer profile can be reconstructed when needed. It should be noted 
that the value of is not changed when the 77 -scale is expanded. This is 
necessary to avoid numerical instability. Using the variable 77 -scale, only 
about 75 points are needed in the 77 -direction. 


Separation is taken to occur at the point at which the wall shear vanishes. 
This point must be located by extrapolation rather than interpolation, because 
the iteration usually fails to converge downstream of the separation point. 


The computation is initiated at each instant using the Hiemenz stagnation- 
point profile (Ref. 12). The turbulent boundary layer downstream of reattach- 
ment of a leading-edge bubble is started from an equilibrium profile with a 
small but finite shear, taken from Ref. 13, matching the computed displacement 
thickness with that of the equilibrium profile. 


In the analysis of the laminar boundary layer, the solution at each Sm is 
used to determine whether transition occurs, accounting for the effects of free- 
stream turbulence and pressure gradient. The following formula, which is a 
modification of a result given in Ref. 14, is used to compute transition 
Reynolds number Re* : 

u tr 


3.6 




fp(A) Re^ 


9860 = 0 


(16) 


where u' is the root-mean- square fluctuation of the free stream due to turbu- 
lence, A is the Karman-Pohlhausen pressure-gradient parameter, A = - ( S 2 <9 p/d x)/ 

fi q e , and Re* is the maximum Reynolds number, based on boundary layer thick- 
er 

ness and local external flow, for which the flow remains laminar. The function 
fp (A) 9 plotted in Figure 4, was derived from data given in Ref. 12 pertaining 
to the effects of pressure gradient on transition, starting from the plot of 
Reg. vs. A (Figure 17.3 of Ref. 12), where Regf is Reynolds number based on 

displacement thickness at the point of instability. Using the Karman-Pohlhausen 
integral method and the plot of Re^ tr -Re^vs. mean pressure gradient parameter 
K (Figure 17.9 of Ref. 12), where Re* and Re^. are Reynolds numbers based on 

momentum thickness 6 at the points of transition and instability, respectively, 
and K is the integral average of - (0 2 <9 p/dx)/^ q e , Reg was computed as a func- 
tion of A and the curve of Figure 4 was constructed. The terms of Eq. (16) 
which account for free-stream turbulence (i.e., Eq. (16) with f p = 1) were 
taken directly from Ref. 14. The form of the terms was derived in Ref. 14 
from arguments concerning the relation of viscous shear to transition. The 
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relation agrees well with measured transition Reynolds number for a flat plate. 
Similar arguments were used in Ref. 14 to derive the effects of pressure 
gradient on transition. The result is an equation of the same form as Eq. (16), 
but with fp simply varying linearly with A. It was felt that the relation so 
obtained was not suitable for the problem at hand, in that it did not give 
good correlation with experiment and does not properly reflect the stabilizing 
influence of a favorable pressure gradient. The term accounting for effects 
of pressure gradient was therefore derived by the procedure previously 
described. 


Viscous Mixing and Reattachment Regions 

For three of the commonly observed types 'of separated flow on airfoils, 
namely short (leading-edge) separation bubbles, leading edge stall (burst 
short bubble) and trailing edge stall experiments have demonstrated that the 
separated shear layer is turbulent at the point where the layer begins recom- 
pressing near reattachment. In the latter case the entire separated flow is 
turbulent because transition occurs upstream of the separation point. In the 
first two cases the shear layer is laminar for a short distance downstream of 
separation, but experiments seem to indicate that transition occurs over a 
distance which is much smaller than either the length of the constant pressure 
mixing region or the length of the fully turbulent reattaching shear layer. 

In these cases then, the assumption of a discontinuous change from laminar to 
turbulent flow in the shear layer seems justified. In the case of thin air- 
foil stall (long bubbles), on the other hand, the transition process in the 
shear layer appears to require a length scale which is not small compared to 
the bubble length. In this case reattachment may begin with the layer still 
in a transitional state. A meaningful flow model for reattachment of a transi- 
tional shear layer would be exceedingly difficult to construct, at least with 
presently available experimental data, and, therefore, thin-airfoil stall will 
not be considered. 

Reattachment of a turbulent shear layer .— In the reattachment region in 
which the bubble thickness shrinks to zero and the streamwise pressure gradient 
is positive, the separated shear layer may reattach onto the airfoil (short 
bubble) or attach onto the layer shed from the lower surface at a stagnation 
point in the wake (leading-edge and trailing-edge stalls). In these flows it 
is assumed that the separated turbulent flow has a wake-like behavior and that 
the equilibrium wall layer usually present in attached turbulent boundary 
layers can be ignored. The wall shear stress is also assumed negligible com- 
pared with the properly normalized rate of change of momentum thickness of 
the layer and the streamwise pressure gradient. These assumptions and the 
analytical model which will be employed here have been used by Todisco and 
Reeves (Ref. 15) and Hunter and Reeves (Ref. 16) to treat supersonic separated 
and reattaching turbulent flows. 
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For negligible wall stress and zero lateral pressure gradient normal to 
the dividing streamline of the bubble the momentum integral and first moment 
of momentum equations are (Ref. 15) 



+ 



8 * 

+ (2H+1) — 


du e 

dx 


d 8 * * dj dH 3J 8 * du e 

^ dx dH dx U e dx 


~ 0 
Kg HR 


Because J = J(H),R = R(H), K^=K^(H) are known profile functions obtained 
from a one-parameter family of equilibrium or self-similar solutions the de- 
pendent variables are the velocity profile shape parameter H = 6/8*, displace- 
ment thickness 8* and the local inviscid velocity at the edge of the turbulent 
layer u e . For supersonic separation and reattachment a third differential 
equation involving these three dependent variables is obtained by coupling the 
integral continuity equation with the Prandtl-Meyer relationship and this third 
equation completes the set. In the case of subsonic separation bubbles, the 
interaction of the separated flow with the outer inviscid stream is much more 
complicated because of the absence of a local relationship between turning angle 
and pressure. In general, a complicated solution procedure is required in which 
the entire bubble shape must be iterated upon until the mixing rate in the shear 
layer and the subsequent rate of pressure rise at reattachment are compatible 
with the induced inviscid flow over the bubble. Such a point by point matching 
is beyond the scope of this study. Instead, approximate relationships between 
certain quantities appearing in the integral equations are obtained from the 
well developed theory of supersonic interactions, and the matching is performed 
in terms of only a few parameters. 


The momentum integral and first moment can be rearranged and combined to 
give the following explicit relations for the derivatives of displacement thick- 
ness and inviscid velocity 
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These equations are solved using the method of successive approximations. 
By first assuming that reattachment occurs over a sufficiently short stream- 
wise distance (compared with the local displacement thickness), so the second 
term on the RHS of the first of these equations is negligible, one obtains 
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s*/s b * = 



f (H) dH 


where 


f(H) = 


-3J + (2H + 1) dJ/dH 


J (H- 1) 


and subscript b denotes conditions at the beginning of the reattachment pres- 
sure rise. Substituting this first approximation for 8* into the differential 
equations for 5* and u e and integrating gives 
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where subscript denotes conditions at the reattachment point (or wake stag- 
nation point if bubble closure occurs in the wake). It is now assumed that 
the streamwise rate of change of the turbulent shear layer velocity profile 
(or rather the inverse, d(x/SR*)/dH i s a universal function for all incompres- 
sible turbulent reattachment processes. This function was evaluated from 
several solutions of supersonic reattaching shear layers with the results shown 
in Figure 5. The variables x and 8 [* are the "stretched" values after having 
been transformed by a modified Stewartson compressibility transformation. 

Three curves for (x-xr)/£* R versus H/Hr for Mach numbers varying between 1.3 
and 2.2 are given and it is evident that the curves collapse into a single 
curve as the Mach number is reduced. Thus, it is assumed -that the curve shown 
for M e = 1.3 defines the streamwise rate of change of the profile shape in all 
turbulent reattachment processes. Values of the derivative d(x/£>R* )/d(H/HR) 
obtained from this curve and values of the profile function obtained from 
equilibrium turbulent boundary layer solutions are given in Table I. 

TABLE I 

PROFILE FUNCTIONS FOR TURBULENT SEPARATED FLOWS 3 


h / h r 

/ x— x r \ / H \ 
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.890 
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1.18 
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NOTEi Hr — Hg — .429, Jj^ — Jg — .654, Rr — .463, and — 

.05 - .0132 (h/Hr). Also, it is^convenient to use the 
product RH as H — * 0 rather than R because in this limit 
R— » oo but RH-»0 . 
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Figure 5 also shows the variation of (x b - x R ) / 8 * , b with H/Hr, that is , 
the reattachment length normalized by displacement thickness at the beginning 
of reattachment (only the curve for M = 1.3 is shown) . This curve can be rep- 
resented quite well by the parabola: 

(X R - x b )/S b * = 10.5 [ 1 - (H b /H R ) 2 ] (19) 


Because and 8^* are determined by the solution of the constant pressure mix- 
ing region this relationship can be used to find the length of the reattachment 
region without having to first compute the velocity and displacement thickness 
distributions through the reattachment region. 


The turbulent mixing region. — For the region between separation and the 
beginning of reattachment in the case of trailing-edge stall, and for the re- 
gion between transition and the beginning of reattachment in the case of 
leading-edge stall, it is assumed that the pressure is constant and that the 
momentum integral and first moment equations reduce to the following: 
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Now, let the subscript t denote conditions at the transition point in the shear 
layer. If the layer is turbulent at the separation point then H t = H s , 8 Z * = 
a s * etc. 

The solution of the first of these equations is 


8 * H = 8 * H t = constant 


and substitution of the first equation into the second gives 
H I dH H j dx * 

Rearranging this expression and integrating from the transition point (or 
separation point for trailing-edge stall) to the beginning of reattachment 
gives 
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and from the solution — constant 


V (H t /H s> 

= OV hJ 

At the transition point \*;e assume that the velocity along the dividing stream- 
line of the separation bubble and the momentum thickness are continuous so 

that* 

^t^lam “ ^Pturb 


and 



Consequently, since H s 6 / 8 * 
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Because x t , (H t /H s ) Iam and (S t * )lant are known from the solution of the 
laminar shear layer, the above expressions are sufficient to determine the 
length of turbulent mixing region and the displacement thickness at the be- 
ginning of reattachment in terms of the unknown profile shape parameter Hb 
at this point. 

The laminar mixing region . — For the laminar shear layer extending from 
the separation point to the transition "point" inside the separation bubble 
the momentum integral and first moment equations are (for zero pressure 
gradient) 



Solution of this pair of equations has been given in Ref. 17. The solution 
presented in that paper was for the variation of the displacement thickness 
and dividing streamline velocity as functions of the non-dimensional length 
scale downstream of separation *2 / 02 R e ^, Since at a laminar separation 

point Hs = 0 s /S g * = .25 the above length scale can be cast into the form 


This procedure for finding the turbulent momentum thickness, displacement thickness and profile shape downstream 
of transition is analogous to the method used by Truckenbrodt for attached boundary layers (Ref. 12). 
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Re x /(Reg * ) 2 = .062 x 2 / 6 S 2 Re x 


Also, H/H s is a function of u* (the ratio of the velocity along the dividing 
streamline to u e ) , so that the solution given in Ref. 17 can be used here with 
only slight modification. The variation of H/H s and S*/5 s * with distance down- 
stream of separation is shown in Figure 6. Here x is the distance downstream 
of the separation point. 

Tentatively, we assume that the length of the laminar shear layer is 
determined by Von Doenhoff's criterion, namely that (R e __ ) = 50 000.* Thus, 

x t 

(JL/H 0 ) and (5*/S *) are known functions of Reo *, as indicated in 
c s lam. c s lam. °s 

Figure 6. The characteristic velocity in each of these Reynolds numbers is 
the local inviscid velocity at the separation point. 

Pressure calculation procedure . — For this part of the analysis, the loca- 
tion of the separation and reattachment points and the pressure at reattachment 
are specified, and it is required to compute the distribution of pressure be- 
tween x s andx R . The basic problem is to determine and x^, • With known, 
the rise in pressure from x^ to x R can be computed directly from Eq. (18). 

Since the pressure at xr is specified and the pressure is constant from x g to 
x b , the pressure distribution is then completely determined. 

The two relations used to determine x b and H b are Eqs. (19) and (20). 
Equation (20) relates x^- x c (or x b -x s , in the case of trailing-edge stall) to 
Hfe, while Eq. (19) relates x R -x b to H b . A simple iterative calculation is 
performed whereby successive values of are assumed and x^, is computed from 
Eqs. (19) and (20). A solution is obtained when the two results for x b agree 
within a prescribed amount. 


Leading-Edge Bubble 

In the analysis of the leading-edge bubble, it is assumed that the begin- 
ning of the reattachment region and the point of transition are coincident. 

The formulations of the viscous mixing and reattachment regions can then be 
employed , as follows . 

From the curves of Figure 6, the conditions at transition, and hence at 
x b , can be obtained, given 8* and 6 at separation from the boundary- layer 
analysis. Since x^ , which is assumed equal to x t is known from Von Doenhoff’s 
criterion ((x t ~x s )q e f v = 50 000), x^ can be computed directly from Eq. (19). 


Gault (Ref. 18) has performed a large number of experiments in which (Re^.) ranged between 25 000 and 75 000 for 

over 80 percent of the data. Consequently, as Gault points out, a unique value for this transition Reynolds number 
is only a rough approximation and is influenced by the turbulence level, pressure gradient, etc. 
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Also, the pressure rise possible in the reattachment region can be obtained 
from Eq. (18). 

It is assumed that the pressure rise across the bubble, from x s to x R , is 
not affected by the presence of the bubble. Given x s and the potential- 
flow solution is then used to compute that increase in pressure. If the pres- 
sure rise computed from Eq. (18) exceeds the required pressure rise computed 
from the potential-flow solution, it is assumed that the mixing process 
accommodates to the lower required pressure rise and that the flow reattaches. 
If the required pressure rise exceeds the pressure rise computed from Eq. (18), 
it is assumed the leading-edge bubble has burst and the airfoil is undergoing 
leading-edge stall. 
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COMPUTATION PROCEDURE 


The procedure used in analyzing the interactions of the flow elements is 
indicated schematically in Figure 7 . The computations are carried out by 
forward integration in time, starting from steady, attached flow. At each 
time step, when the flow is attached, prescribed airfoil motions are used to 
derive the potential flow. The flow external to the boundary layer is then 
computed. The boundary layer and leading-edge bubble are then analyzed to 
determine whether the leading-edge bubble has burst (leading-edge stall) or 
the turbulent boundary layer has separated (trailing-edge stall) . If the flow 
remains attached, time is incremented, the foil motions are again prescribed 
and the procedure repeated. 

If the airfoil stalls, the location of the separation point is assumed, 
the point of reattachment and pressure at reattachment are computed by a 
method to be described subsequently, the pressure in the dead-air region is 
computed and the potential-flow solution is derived from the prescribed air- 
foil motions and dead-air pressure distribution. The boundary layer is then 
analyzed and the computed separation point location is compared with the 
assumed location. If the assumed location differs by more than a prescribed 
amount from the computed one, a new estimate of the separation point location 
is made and the potential flow and boundary layer are again analyzed. When 
the assumed and computed separation point locations are in satisfactory 
agreement or when four iterations have been performed, time is incremented. 

The length of the dead-air region when the airfoil is stalled is computed 
by deriving a bubble growth rate from the potential solution and integrating 
in time. The growth rate is computed as follows. The rate of increase of 
mass of the dead-air region is given, to first order, by 



x s 


The total mass of the_dead-air region is roughly proportional to pl s Y max , 
where 1 = -r and is the maximum of Y(x) , 

s ix. s max 
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The primary justification for computing x R in this way is that it gives a 
closed bubble in the limit of steady separated flow. At the onset of stall, 
when there is no potential-flow solution from which to derive a growth rate, 
l s is assumed to equal U . 

When the dead-air region terminates doxmstream of the trailing edge, the 
pressure at x R is taken to be the free-stream pressure. If x R is less than b, 
as occurs at the onset of leading-edge stall, the pressure at x R is taken to 
be the pressure at that point if the airfoil were not stalled. 

When the airfoil has undergone leading-edge stall, there must be some 
means of determining when the airfoil regains attached flow. Closely con- 
nected with this determination is the process by which reattachment occurs. 

It was assumed that when an airfoil regains attached flow after having under- 
gone leading-edge stall, the dead-air region is washed off the airfoil from 
leading edge to trailing edge at the free-stream speed. To determine whether 
this process should be initiated, the following procedure was first attempted. 
It was postulated that during the previous time step, the airfoil had unstalled 
and the dead-air region was partially washed off the airfoil. After computing 
the potential flow with x s so positioned, the boundary layer and leading-edge 
bubble were analyzed. If the leading-edge bubble burst, the postulated unstall 
could not have occurred and the analysis would proceed with the separation 
point restored to a position near the leading-edge. Unfortunately, this pro- 
cedure caused difficulty with the boundary-layer analysis, because at large 
angles of attack the displaced dead-air region generated a local area of un- 
favorable pressure gradient near the stangation point, probably because of 
the inaccuracy of the linear potential-flow model used. The procedure was 
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therefore altered somewhat by removing the dead-air region entirely from the 
flow prior to analyzing the boundary layer and the leading-edge bubble. Of 
course, if it is found that the leading-edge bubble does not burst and unstall 
should be initiated, the washing off of the dead-air region is still carried 
out. 


RESULTS OF COMPUTATIONS 


Calculations have been carried out of the loading on an airfoil during 
transient and sinusoidal pitching motions and during unsteady stall induced 
by a series of discrete vortices convected past the airfoil at the free-stream 
speed. Two different airfoil sections were analyzed. These sections, which 
are sketched with their offsets in Figure 8, were the ones used in the tests 
reported in Ref. 3. The symmetric section is an NACA 0012 section modified 
by a small trailing-edge extension. The cambered section, designated V-23010 
in Ref. 3, was formed from a 0012 section by cambering the leading-edge portion. 

All calculations were performed for a chordal Reynolds number of two mil- 
lion (Re^ = 10^). The measured static variation of lift and moment, from 
Ref. 19, indicate that the sections analyzed undergo leading-edge stall at 
that Reynolds number. Preliminary calculations of the loading during transi- 
ent pitching through stall confirmed this. Therefore, in order to conserve 
computer running time, the occurrence of trailing-edge stall was precluded 
by omitting the analysis of the turbulent boundary layer downstream of the 
leading-edge bubble. Thus, the possibility of unsteady effects suppressing 
leading-edge stall to the extent that trailing-edge stall dominates were not 
taken into account in the computations. However, there is no indication in 
the data of Ref. 19 that trailing-edge stall occurs under the conditions 
analyzed. 


Transient Pitching Motion 
and Static Stall Characteristics 

Computations of the loading resulting from transient pitching motions 
have been carried out for both airfoil sections. Pitch angle was increased 
linearly with time up to a selected value and then held constant until steady 
loading on the airfoil was established. The variation of pitch angle, normal- 
force coefficient, moment coefficient and length of the dead-air region with 
time, for a typical case, are shown in Figure 9. 

Of particular note in Figure 9 is that the instantaneous normal force 
exceeds the maximum static value computed by a considerable margin at stall. 
There are two effects contributing to the overshoot. One of these is a redis- 
tribution of loading over the chord, caused by the pitch rate, which reduces 
the suction peak and pressure gradient near the leading edge. This can be 
seen from a comparison of various flow quantities for static and dynamic 
loading below stall as listed in Table II. The quantities for unsteady 
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TABLE I! 


COMPARISON OF STATIC AND DYNAMIC LOADING BELOW STALL 


Quantity 

Static 

Dynamic 
(b<?p/u= .025) 

6 p -deg. 

11.5 

12.9 

C n 

1.252 

1.354 

q e /u 

max 

2.985 

2.892 

q e /u 

S 

2.832 

2.742 

Length of leading 
edge bubble-semichords 

.0256 

.0265 

Increase in 
Cp possible 

2.095 

1.964 

Increase in 
Cp required 

2.093 

1.949 


loading were taken from the run plotted in Figure 9 at the instant when 
0p = 12.9 deg. The normal force in the unsteady case is seen to exceed the 
static normal force, even though the maximum flow at the surface, and hence 
the suction peak, is considerably less for the dynamic than for the static 
loading. As can be seen from the comparison of the pressure distributions 
for the two cases shown in Figure 10, the additional normal force is gener- 
ated over the aft 90 percent of the chord on the pitching airfoil, the load- 
ing near the leading edge being slightly less than on the static airfoil. 

This is caused by the pressure distribution associated with the pitch rate, 
the boundary condition for which is equivalent to parabolic camber on an air- 
foil in steady flow. The effective camber shifts the center of pressure rear- 
ward, reducing the suction peak. As can be seen from Table II, the reduced 
flow magnitude at separation in the dynamic case makes the separation bubble 
longer and decreases the amount of pressure rise in the bubble which is possi- 
ble, but this is more than offset by the reduction in pressure gradient, which 
lowers the required pressure rise even more. The reduction in pressure grad- 
ient in unsteady flow was in fact postulated by Carta in Ref. 20 as a mechan- 
ism for dynamic overshoot, without considering the direct effect on the leading- 
edge bubble. 

The other effect contributing to the increased normal force derives from 
the flow which results when the dead-air region terminates on the airfoil. 

When the reattachment point is in the vicinity of midchord, the presence of 
the bubble induces a substantial increase in loading over the aft portion of 
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the airfoil which more than compensates for any loss in lift over the dead-air 
region, resulting in a net increase in normal force. This can be seen from 
the variation of pressure coefficient over the chord plotted in Figure 11 
for three different pitch angles, taken from the run plotted in Figure 9. For 
0 p - 14.3 deg., the dead-air region extends from very close to the leading edge 
to approximately midchord. The loading is seen to exceed that for 0 p =12.9 
deg (just prior to bursting) over the whole aft portion of the airfoil, result- 
ing in a normal-force coefficient of 1.59, compared to 1.35 at 0 p = 12.9 deg. 
When the reattachment point reaches the trailing edge and beyond, the loading 
rapidly drops to near the static variation during stall, as can be seen from 
the variation of Cp y for 0 p = 17.2 deg in Figure 11. 

Of further note in Figure 11 is the similarity of the pressure distribution 
for 6 p = 14.3 deg to that which would be induced by a strong vortex. Vorticity 
which is shed from the leading edge has been postulated (Ref. 4) as a primary 
potential-flow element of dynamic stall. Shed vorticity emanating from the 
leading edge is not a part of the model employed, per se , so a different 
explanation is required. From a mathematical standpoint, the cause of the 
sharp rise in suction downstream of reattachment is the square-root singularity* 
in y required there to satisfy the f low-tangency condition. A physical inter- 
pretation could be that fluid, upon passing from a region of nearly constant 
pressure to one of attached flow, must accelerate in a manner analogous to 
passage over a leading edge, causing a sharp suction peak. 

The variation of static normal-force and moment coefficients with angle 
of attack were developed from the transient pitch computations by parametrically 
varying the maximum pitch angle. The curves for the symmetric and cambered sec- 
tions, and the corresponding measured coefficients from Ref. 19, are shown in 
Figures 12 and 13. The computed and measured static coefficients are seen to 
agree quite well, qualitatively. That the method so accurately predicts the 
maximum normal force on the symmetric section was fortuitous, in light of the 
approximations made in the calculations, and should not be regarded as a 
measure of the quantitative capabilities of the method. The differences be- 
tween measured and calculated coefficients for the cambered section (Figure 13) 
can be attributed in part to a larger camber effect in the computed result, 
there being a considerably larger predicted normal force at a = 0 than was 
measured. Also, the center of pressure is quite a bit aft of the quarter-chord 
in the predicted load, as reflected in the moment variation. Apparently the 
linearized model introduces a relatively large error when used to represent 
camber of this type. 


Sinusoidal Pitch Oscillations 

Computations were carried out of the load on both symmetric and cambered 
sections resulting from sinusoidal pitching about the quarter-chord. Because 
the calculations were initiated from steady flow, it was necessary to carry 
out the computations until periodic loading was established. It was found that 


The singularity results from the linearized boundary condition. A uniformly valid solution v/hich is finite at the 
reattachment point could be derived usin£ Lighthill's approach, but it v/as not needed, because the boundary layer 
is not analyzed in the vicinity of that point. 
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Figure 13 STATIC NORMAL-FORCE AND MOMENT COEFFICIENTS FOR CAMBERED AIRFOIL 





just two cycles of oscillation were sufficient to obtain periodicity. The time 
increment was selected to give 36 time-steps per cycle of oscillation. A typi- 
cal run required about 40 minutes of central processor time, using an IBM 
360/75 computer. 

Symmetric section . — The results of the computations for the symmetric 
section are shown in Figures 14 through 19 as plots of normal-force and moment 
coefficient vs. pitch angle. The static variation of those coefficients is 
shown as dashed curves on the figures. The loading for a reduced frequency 
k (k = co b/U, where co is circular frequency) of .13 is given in Figures 14, 15, 
and 16 and k = .26 for Figures 17, 18 and 19. The amplitude of oscillation 
for all cases is 4.5 deg. The three mean pitch angles for each reduced fre- 
quency are 9, 13.5, and 15.5 degrees, as indicated on the figures. 

For purposes of comparison, the measured loadings shown plotted in Fig- 
ures 20 through 25 we re obtained from the data given in Ref. 19. The param- 
eters for those plots correspond nominally to those of Figures 14 through 19, 
respectively. The measured coefficients were computed from the harmonic 
analyses given in Ref. 19 of data taken from tests at a free-stream Mach 
number of .2. 

It is noted, first, that the measured loads for k = .13 (Figures 20, 21, 
and 22) exhibit considerable dynamic overshoot of C n , while none is evident 
in the computed results at that frequency (Figures 14, 15, and 16). The com- 
puted moment coefficient is in good agreement, qualitatively, with the measured 
coefficient, however. The unstable moment variation (loops traversed clockwise 
in a plot of C m vs. ) are in clear evidence, particularly at a mean pitch 
angle of 9 deg (Figure 14). Further, the variation of normal force during the 
process of flow reattachment exhibited in the measured loading is well repro- 
duced in the computed result, there being an initial sharp rise in C n as the 
dead-air region is being washed off, followed by a slight drop in C n after the 
airfoil regains attached flow. 

At the higher reduced frequency, the computed moment coefficients (Fig- 
ures 17, 18, and 19) continue to exhibit the unstable variation in evidence 
in the experimental result (Figures 23, 24, and 25). Also, the computed 
normal-force coefficients for mean pitch angles of 9 and 13.5 degrees (Fig- 
ures 17 and 18) are seen to exceed the maximum static value by a considerable 
margin. It would appear that overshoot is only detected at the higher fre- 
quency at least partly because of the substantial errors introduced by the 
linear model in computing the flow near the leading edge, because it was found 
in the analysis of transient pitching that bubble bursting is very sensitive 
to changes in the flow there. 

In further comparison of the results at k = .26, it is noted that the 
measured C n does not decrease as rapidly as the computed load at the onset of 
stall, making the experimental C n loops considerably narrower than the computed 
ones. Because the decrease in computed normal force does not occur until the 
reattachment point has progressed to the vicinity of the trailing edge, it 
would appear that the rate of growth of the dead-air region, assumed to be 
equal to the free-stream speed in the analysis, is in fact a good deal less 
than the free-stream speed. 
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For the case computed with the highest pitch angle at k = .26 (Figure 19), 
it can be seen that the airfoil remained stalled throughout the cycle, while 
the measured results at the same frequency and comparable pitch angle (Fig- 
ure 25) indicate the airfoil did unstall. It is interesting to compare 
Figure 19 with the measured load for a case where the airfoil did remain 
stalled. The result for such a case, also with k = .26, is plotted in Fig- 
ure 26. Note in Figure 26 that the airfoil developed a substantial increase 
in load while pitching up, but the computed results show only a slight in- 
crease in C n over that part of the cycle. This indicates that there is a 
significant variation of pressure in the viscous mixing region caused by 
unsteady effects which are not being taken into account in the analysis be- 
cause a quasi-steady model of that region was used. 

Cambered section . — It was found in Ref. 3 that the cambered section gen- 
erally exhibited smaller excursions (less negative values) of C m and reattach- 
ment started earlier in the cycle, compared to the results for the symmetric 
section. The load on the cambered section was computed to determine whether 
these effects could be detected. The results for reduced frequencies of .13 
and .26 are plotted in Figures 27 and 28, respectively. Both cases have a 
mean pitch angle of 15.5 degrees and pitch amplitude of 4.5 degrees. Because 
the stall angle of attack of the cambered section is considerably larger than 
that of the symmetric section, Figures 27 and 28 are most closely comparable 
to Figures 15 and 18, respectively. The moment coefficients of the cambered 
section appear to be about the same magnitude as those of the symmetric sec- 
tion. However, reattachment for the cambered section does occur significantly 
earlier in the cycle, in agreement with the findings of Ref. 3. Note that 
substantial overshoot of C n again occurs only at the higher frequency. 


Wake-Induced Stall 

The method was applied to the analysis of wake- induced dynamic stall of 
a rotor blade during a maneuver. Large high-frequency oscillations in dif- 
ferential pressure Ap, torsional moment TM, lift and aerodynamic moment on the 
retreating side, with blade azimuth angle ^between 270 and 360 degrees, were 
detected in flight test data, as shorn in plots of Figures 29 and 30, taken 
from Ref. 21. It was asserted in Ref. 21 that this response was the result of 
dynamic stall induced by previously formed tip vortices which, under the 
maneuver flight condition, pass under the blade at the azimuth positions indi- 
cated in Figures 29 and 30. To test this hypothesis, calculations were carried 
out of loading resulting from passage of a two-dimensional airfoil over a series 
of three equally spaced discrete vortices. The airfoil was free to respond in 
pitch about its quarter-chord point, with spring restraint and inertia chosen 
to give a reduced natural frequency k# = ogb/U of .382, which corresponds to a 
pitch natural frequency of about 40 Hz for the parameters of the flight test . 
Pitch angle was computed from 
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Figure 29 MEASURED DIFFERENTIAL PRESSURE TIME HISTORIES FOR THE 95% 
BLADE RADIUS IN A 1.5g PULLUP MANEUVER AT AN ADVANCE RATIO OF 0.22 

(from Ref. 21) 










0p(t) = d Q + [0 p (t- At) - e o ] COS ((00 At) + 


sin {(£>q At) 


flp(t-At) 

co e 


2 mj 

+ — [2C m (t — At) - C m (t— 2At)] [l — cos (tog At) ] 

k 0 2 
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k 0 3 U c At 


0p (t) = (t- At) cos (coq At) - <oq [0p (t - At) - 6 q ] sin (coq At) 

2 mj 


+ COQ c m (t - At) sin (coq At) 


4 

where 0 Q is the angle of zero restoring moment, and mj = pb /Iq , Iq being pitch 
moment of inertia per unit span. A value of .03583 was used for mj in the cal- 
culations. The above formulas were derived by integrating the pitch_equation 
using a linear extrapolation for c m . Dimensionless vortex strength r = T t /ub 
was assigned a value of .8636, which corresponds to the strength r t of the tip 
vortex formed by a heavily loaded blade. The airfoil section of the rotor 
blades used in the flight test is a 0012 with a trailing-edge extension, so 
the symmetric section employed in the previous calculations was used for these 
as well. 

A total of nine cases were run, parametrically varying streamwise and 
vertical spacing of the three vortices to determine whether one or more com- 
binations produce load variations similar to those of Figures 29 and 30. The 
variations with time of pitch angle (comparable to torsional moment in Fig- 
ure 30), lift coefficient and moment coefficient for those cases are shown 
plotted in Figures 31 through 39. 


— [C m (c- At) - C m (t - 2 At)] [l - cos ( coq At) ] 
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It should first be noted, for purposes of comparison, that the computed 
oscillations in pitch angle are at approximately the pitch natural frequency 
(18 semichords traveled per cycle) regardless of the vortex spacing, and the 
measured oscillations in torsional moment, which correspond directly to blade 
pitch response, are also at about the pitch natural frequency. It. can be 
seen that, while the three cases with vertical spacing of four semichords 
(Figures 37, 38, and 39) exhibit considerably less loading variation than the 
measured results, the computed loading for the three cases with vertical spac- 
ing of one semichord (Figures 31, 32, and 33) as well as those with two-semi- 
chord vertical spacing (Figures 34, 35, and 36) are quite' similar to the meas- 
ured variation of Ci and C m in the interval of interest. Lift coefficient 
varies from about .6 to 1.3 and C m varies from near zero to about -.15 in both 
computed and measured results. The wave forms are also similar, particularly 
for the largest streamwise spacing of 16 semichords (Figures 33 and 36). The 
rapid variations in loading which occur during reattachment are absent from 
the measured loads, but it is doubtful that the instrumentation employed for 
the flight test would have detected oscillations of such high frequency (about 
six times the pitch frequency, or 240 Hz). 

The computed differential pressure coefficients at selected chordwise 
stations for the runs with streamwise spacing of 16 semichords and vertical 
spacings of one and two semichords are plotted as a function of time in 
Figures 40 and 41, respectively. Comparing these figures with Figure 29, 
shows that the variations in computed chordwise loading are also similar to 
the measured pressure variations. In particular, the computed loading repro- 
duces the phase lag, or time delay, of the loading at the aft chordwise 
stations with respect to the forward ones. 

It can be concluded from the results obtained, then, that the wake- vortex 
mechanism hypothesized in Ref. 21 is causing the large high-frequency response 
of the blade. The results further explain why the oscillation is at the pitch 
natural frequency, rather than at the excitation frequency corresponding to 
the vortex spacing. The large nose-down moment exerted on the airfoil when 
it stalls causes the blade to rapidly pitch down and unstall. It was found 
in the computations that the blade only travels one or two semichords between 
stall onset and the beginning of flow reattachment. Thus, the blade is effec- 
tively excited by a series of discrete impulses, causing it to t oscillate at 
its natural frequency. 
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CONCLUSIONS 


The method developed is capable of reproducing the essential features of 
dynamic stall, as demonstrated by the good qualitative agreement obtained be- 
tween computed and measured loads on a pitching airfoil. The results indicate 
that dynamic overshoot of the normal force is caused by a redistribution of 
pressure associated with the pitch rate as well as loading induced on the aft 
portion of the airfoil by the presence of a dead-air region on the forward 
portion. Quantitative differences between theory and experiment can be attrib- 
uted partly to the use of a linearized representation of the potential flow. 

The assumption that rate of growth of the dead-air region at the onset of 
leading-edge stall is the free-stream speed and the use of a quasi-steady 
model for the viscous mixing region also cause quantitative differences between 
the computed and measured loading. 

The large torsional response of helicopter blades during a maneuver which 
had been detected in flight tests can be attributed to dynamic stall induced 
by previously formed tip vortices, as hypothesized in Ref. 21. The response 
is primarily at the torsional natural frequency, rather than at the excitation 
frequency, because the severe nose-down moments in stall cause the blade to 
rapidly unstall, effectively loading the blade with a series of impulses. 
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APPENDIX A 


SOLUTION FOR STALLED FLAT PLATE IN STEADY FLOW 


First, consider the case with x R > b . Equations (2) through (5) become, 

for steady flow with constant pressure p in the dead-air region: 

d 


*£ 


y<0 

X- $ 


= 2 U a, — b < x < x 


(A-l) 


a(x) + 

77 


£ 


i I yUU£ 

x- £ 


= 2U a, x < x < b 


(A-2) 


y(x) + 

7 T 


{ 




a(£)d£ _ 2 <Poc-Pd) 

X - ^ pu 


Xg < X < b 


(A-3) 


^ f. 




r(() d£ 2 <Peo-Pd) 




pU 


b < x < x- 


: R 


(A-4) 


The first step is to reduce the number of equations by solving for y on 
the interval (-b, x s ) in Eq. (A-l) and solving for a on the intervals, x R ) 
from Eq. (A-4). Formally inverting those two singular integral equations (see 
Ref. 22), it is found that 
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where A' and B' are as yet undetermined constants. If Eqs. (A-5) and (A-6) 
are substituted in Eqs. (A-2) and (A-3) , respectively, and certain of the 
resulting integrals are evaluated, the following pair of integral equations 
is obtained: 





(A-8) 
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It is required that o be well behaved at x = x g , which, from Eq. (A-7) , allows 
the value of A' to be assigned: 


A ' 


= 2 U a x s + 




y(jj) d 77 


Equations (A-7) and (A-8) can be combined by formally solving Eq. (A-8) 
for a : 
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If a is to be at least integrable at x=b , the quantity in brackets must vanish 
as x-» b . This provides a relation between the undetermined constants C and B 
The continuity of a at x = x s provides the other relation needed to solve for 
those constants, with the result that 
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V x ~* s 

Xr-X 
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Substituting Eq. (A-9) and the relation for A' in Eq. (A-7), it is found that 
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Now, let 




x + b 

g GO = (x R -x) ^/- — — y(x) 


z(x) = 




x + b 


With some manipulation of the integrand of Eq. (A-10) , that relation can be 
written in the form 
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This equation can be solved for g . Letting z g = z (x s ) and z^ = z (b) 
Eq. (A-ll) gives 
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where the undetermined constant was assigned to make g well behaved at z^ and 
k a = U a, kp = (pc* - Pd)/pU . If the integrand of Eq. (A-12) is expanded in partial 
fractions, the resulting integrals can be evaluated by standard methods. By 
substituting the result in previous relations, expressions for y and <? over the 
whole interval for which they are defined can be obtained. The complete solution 
for x R > b is: 
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while C Q and have the same definitions as C Q and Cj , respectively, but 

with z replaced by r o and z, replaced by r_ . 
s s b K. 
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APPENDIX B 


POTENTIAL FLOW FORMULATIONS 


The coefficients and inhomogeneous terms of Eqs. (8a) and (8b) are given 
by the follov 7 ing expressions* For simplicity all lengths in this appendix have 
been made dimensionless using semichard as reference length. It is understood 
that any complex term is omitted. 
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[en|l + a| -1] 

- ^ 2 ^ ■ [j n I x - b I - 7] + " ' 2 ' ^ " [ fc I 1 + b I - y \ 

— (b — a) -j( x — b) jjEn |x — b| — lj + (1 + b) j^n | 1 + b | — lj | 


% = f p < x a k ) 


* U o 


4A£— (1 + V 


1 - 3 x., + 


3<1 -x 2) 


«k 24f„ 


A„ (t-At) + 


y Aj(t-At)J 


-A 0 (t-2At) +— Aj(t-2At) 


U„ 


UAt 


2I y (% k , t-At) + 2 I ff (S v t-At) 

- 7 V ( V t_2At) “7 c - 2At) 


for x < 1. 

a k 
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the terms involving A Q and A i and the ly T s being dropped for 


l a (x, t) 

Iy (X , t - 


with x 
For -1 

V x) = 


= B Q (t) ( x i c ) + 



i = 2 


g B . (x,t) B| (t) 


v&t) =(n-Q + sin 6 ) |A 0 (t - v At) + y 



1 

2 



sin (n + 1) 6 
n + 1 


sin (n - 1) 6 
n — 1 


A„ (t - vAt) 


2(1 - x R ) ^ 2 j^(x R - x s ) (1 - x ) ( x - x R )] ^ B 0 (t - i/At) 

- — y w (t) (1 + x) (1 - x 2 ) , 1 / = 0, 1, 2 ; 

4 1 + v 

- COS0 . 

< x <1, with x = cos 6 9 


2 (Poo-Pd ) 

pUU o 


2U 

™o 


max 


C 1 

— + (1 — cos 0) 
2 


Nf 


Z 

n = 1 


t n cosn0 


2UT, 


max 





nsinn 6 (1 — cos 0) 
sin 0 


— cosn 0 


'n 


while for x > 1, 



Nr 


f (x) = 

p puu 0 


2 (P M -P d > 2 UT„ 


UU„ 




n = 1 


Nr 


2 U T_ 


LL 


where I D = 


) Jx (n + 1 ) - n I n - I 


Pn Pn — 1 


n =b 1 


G - v 2 Z )” 

yP~l 


The flow at the surface and the pressure coefficient are given by the 
following expressions (again letting x = cos 8 ) : 


% (x, Oi , t) ± 1 


Uo 

where 


yr+x+ r /2 


o' - L 


— -JlT* + f q (x, Oi, t) 


Nr 


7 T"f q (x, 0 ± , t) = 2 T max cos 0/2 

U o 


1 I— 

E n sin n 
sin 0 


(1 — cos 0 ) - cos n 6 


n = 1 


N. 


t A Q (t) sin 6/2 ± cos 8/2 i 

n = 1 


A n sin n 6 


- 3/2 


+ — (1 + x) ( 3 x - 1 ) (t) + (1 - x R ) 


(x R- x s^l 
. X -*R 


-x) 


1/2 


(l + 3x R - 4x) B 0 (t)J 


N„ 


— cos 6/2 


L Xr>— x 


B 0 (t) + cos 8/2 


E 

n = 2 


fg (x, t) B n (t) 
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Cp (x, 0(i , t) = 2 


q e (x,0-i, t) 

+ 1 

u 


U o (3 1 . ( 

+ - <{— I y (X,t) - 2Iy(X, t- At) + Y Iy(X, t — 2 At)) 


U Z A^ 


U ° 3 1 I 

+ {— Ut)-2 l a (t- At) +— I_(t— 2At) 


u 2 a^ 


O 

U 

u 2 


N f 


+2(1 — cos 0) 


£ 

n = l 


t n cos n 6 


101 




APPENDIX C 


BOUNDARY LAYER FINITE DIFFERENCE RELATIONS 




APPENDIX C 


BOUNDARY LAYER FINITE DIFFERENCE RELATIONS 



K 3 q, 


m + 1 , n 

\ 

d Js 

ds / m + 1 , n 


2 U 0 At J s m + 1 , n 




m+ 1 , n 

m 'Im + l, n ~ b m 1 s mn + c m ^ m _ 1>a 


m ( s m + l _s m) ^ s m + l s m-l^ 
s m + l " S m — 1 


( s m + l- s m) (s m- s m-l> 


s m + l 


(s + i - s m— 1 ) ( s m s m — 1 ^ 


d _% 

d 1 ! /m + 1 , a 


^ s m + l, n + 1 + Cn ^ s m+l, n 


V a - n n -1 




m + 1 , n— 1 


d ” (7 ?n + l ~ ’Jn-l^n+l " > 


n (Vn-Vn-O ^n + l-^fl) 

^n+l ~ % 

fn = tn~ n “ *7,1-1) ^n + l ~ ’Jn-l j 


' d2 % 
dr , 2 


* “n q s m+lj n +2 + qs m+l, n + 1 + Ya K + l, ° 


+ 8 n 


m + 1 , n 


a n = 2 


2 V n ~ Vb+I ~ r 'a+ 1 


( 7 7 n + 2 — ^n — l^^n +2 7 ?n + l^^n + 2 


r ] 


(t — 2 At)l 


qs m + l, n -1 
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2 


»n " 


7n + 2 + 7n + 1 2 7n 


(7 ?n + 2 ~ 7n - 1> <7 n + 1 “ 7„ - 1> (7„ ~ 7n - 1> 




<7n + 1 ~ 7n> : 


(^n - _ l) 3 S n - (7 ?n+ 2 - i? n ) «„ 


Xn _ ( Q n + + 


<?P 


V 5 s /m + 1 \2 U o 2 aJ 


e m +■ 1 


Uo 2 


3 q e , (t) - 4 q e (t - At) + qg ( t — 2 At) 

e r,i +1 e m + 1 m + 1 


£mV„ . - < £ ) + , W ~ % m (t > 

m + 2 m + 1 m 


fn. = 


s m + 1 s m 


( s m + 2 s m ^ ^ s m + 2 s ra + 1 ^ 


(s m + 1 ~ s m > (s m + 2 “ s ra + 1> 


s m + 2 s m + 1 


J'™ = 


^ s m + 1 s m ) ( s m + 2 s m ^ 


The gradient at the wall, needed to evaluate the shear, is given by 


'dq. 


iri +1 ,l- k2 qs m + l,2 k 3 qs m + l,3 + k4qs m + l,4 
73 V4 


k 2 = 


k 3 = 


k 4 = 


7 ?2 (7 /3 _7 72^ 7 ?4 ~ 7 ?2 ) 
72 74 

7 ?3 (, ?4“ , ?3 ) ^ ^ 3 — ^72 > 
72 73 

74(7 4 “ 73) (74 ~ 72) 
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